close all
pointA=[0,0,0];
pointB=[1,1,0];
pointC=[0,1,1];

normal = cross(pointA-pointB, pointA-pointC); %# Calculate plane normal
%# Transform points to x,y,z
x = [pointA(1) pointB(1) pointC(1)];  
y = [pointA(2) pointB(2) pointC(2)];
z = [pointA(3) pointB(3) pointC(3)];

%Find all coefficients of plane equation    
A = normal(1); B = normal(2); C = normal(3);
D = -dot(normal,pointA);
%Decide on a suitable showing range
xLim = [min(x)-1 max(x)+1];
zLim = [min(z)-1 max(z)+1];
[X,Z] = meshgrid(xLim,zLim);
Y = (A * X + C * Z + D)/ (-B);
reOrder = [1 2  4 3];
figure();patch(X(reOrder),Y(reOrder),Z(reOrder),'b');
grid on;
alpha(0.3);
xlabel('x_1');
ylabel('x_2');
zlabel('x_3');

hold on
axis square
rp=[0 1 0];
line([0 rp(1)],[0 rp(2)],[0 rp(3)],'LineWidth',2,'Color',[0 0 1])
angle=(90+30)*pi/180;
t=0.15;
line([rp(1) rp(1)-t*cos(angle)],[rp(2) rp(2)-t*sin(angle)],[0 0],'LineWidth',2,'Color',[0 0 1])
angle=(90-30)*pi/180;
line([rp(1) rp(1)-t*cos(angle)],[rp(2) rp(2)-t*sin(angle)],[0 0],'LineWidth',2,'Color',[0 0 1])
text(-0.1,-0.1,-0.1,'Plane','FontSize',14)
text(0,1.25,0,'r''','FontSize',14)

